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ABSTRACT 



Hydrographic surveys for nautical charting contain many 
discrete data points. Analjrtical models for ocean bottom 
topography could save computer storage and reduce the com- 
plexity of automating the nautical charting process, but 
they must meet stringent accuracy requirements. Polynomials, 
double Pourier series, finite elements, Duchon's analysis, 
Shepard's formula and Hardy's multi quadric analysis were 
investigated as possible modeling techniques. Multiquadric 
analysis in v/hich the surface is represented by an analyti- 
cal summation of mathematical surfaces such as cones and 
h 3 rperboloids v/as the only method found to be suitable. An 
iterative method of model point selection v/as found to give 
the best results. Smooth and unambiguous junctions of 
adjacent models were made by using a Hermite polynomial 
v/eighted sum of overlapping areas. Highl 3 r irregular surfaces 
can be represented by about 20% of the original sur'rey data 
points; more regular bottom topograph;/ can be represented 
by a smaller percentage. 
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I. INTRODUCTION 



A. BACKGROUira 

The ocean bottom is a continuous but generally irregular 
surface. In the deep oceans there are vast areas of abyssal 
plains interrupted by mid-ocean ridges, sea mounts and con- 
tinents, The continental shelves and coastal areas vary 
from smooth flat bottoms to highly irregular surfaces with 
deeply gouged glacial troughs or coral and rock pinnacles. 
Many geological formations which are found on land such as 
canyons, mountains, domes, faults, etc., are also found on 
the continental shelves. The shape of the ocean bottom is 
difficult to determine since it cannot be seen or photo- 
graphed except in very shallow areas and, direct measurement 
requiring occupation of the ocean bottom is costly and often 
impossible. 

There are many reasons for which the shape of the ocean 
bottom must be known. Historically, safety of navigation 
has been the most urgent reason. Nautical charts are com- 
piled from many sources to aid the navigator. These charts 
depict the coastlines and ocean bottom features using con- 
tour lines and selected depths. 

The primary sources of depth data for nautical charts 
are hydrographic surveys. These stirveys represent ocean 
bottom topography by discrete data points which are defined 
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by geographic position and depth below a specified water 
level datura. Until the mid- twentieth century, these depths 
were determined by lowering a weight on a calibrated line 
imtil it touched bottom. The vessel position was usually 
determined by measurements with sextants. Using these 
manual methods, data acquisition was very slow and only a 
minute percentage of the bottom was sampled. There were 
many sources of error in the observational procedures. A 
typical survey had a few hundred data points from which the 
surface shape between points had to be inferred. Data pro- 
cessing was easily handled by manual methods. More recently, 
electronic positioning equipment and depth soimding instru- 
ments have been used in semi -automated and automated systems. 
These systems allow almost continuous sampling of the ocean 
bottom along the vessel track. They have increased the 
accuracy of the data and the completeness of bottom coverage. 
As a result, depths need to be inferred between vessel tracks 
but not along the tracks. A typical survey of this type 
contains between 2,000 and 20,000 data points. These sys- 
tems increased the data acquisition rate to such an extent 
that manual data processing methods could not keep up with 
data acquisition. Computer aided systems for processing 
and verifying the data were developed in the 1960's and 
1970 ‘s. 

Producing a nautical chart requires compilation of 
many hydrographic surveys, shoreline manuscripts, and other 
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dociunents. This remained a manual process until the mid 
1970 's. At this time, the National Ocean Survey (NOS) of 
the United States National Oceanic and Atmospheric Adminis- 
tration (NOAA) began development of a computer assisted 
chart compilation and production system (Moses and Passauer, 
1979 ). This system requires on-line storage and manipu- 
lation of large blocks of discrete point data from hydro- 
graphic surveys. The density of these data from modem 
surveys make this a complex and costly process. 

In an effort to produce one hundred percent bottom 
coverage for critical areas, multi-beam sounding systems 
(Hopkins and Mobley, 1978), airborne ]aser depth measuring 
systems (AYCO Everett Research Laboratory, 1978), and airborne 
water penetrating photography systems (Keller, 1976) have 
been developed. Some of these systems have proved that one 
hundred percent bottom coverage is feasible. They have 
also created another problem concerning representation of 
the data and its use in the compilation of nautical charts. 

The data from the multi-beam sounding systems for a typical 
survey would be equivalent to several hundred thousand dis- 
crete data points. Data from a laser system would be even 
more dense. The photo gramme trie method uses stereographic 
images produced from aerial photographs. This can be con- 
sidered to be truly continuous data, but such data is diffi- 
cult to represent in a digital computer. The usual method 
to represent this data is to select the most representative 
and most critical depths for use as if they were from a 
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conventional siirvey. For a bottom with little relief, 
this method is satisfactory but as bottom relief increases, 
considerable detail and completeness is lost. 

B. MATHMATICAIi MODELS FOR OCEAIT BOTTOM TOPOGRAPHY 

The density of data from modem hydrographic surveys 
has made the automation of chart compilation difficult. A 
possible solution to this problem v/hich is investigated by 
this thesis is the use of a sxirface defined by an analytical 
expression to approximate the ocean bottom topography. 

Such a mathematical model would be used to compute a depth 
at any geographic position within the bounds of the model. 

In order to be useful, such a model must require consi- 
derably less data storage for the parameters which define 
the model than was required by the original set of discrete 
points. 

The accuracy of the model is of utmost importance. The 
United States government can be held liable for vessel 
groundings or accidents at sea which are due to inaccurate 
charts. Special Publication 44 of the International Hydro- 
graphic Bureau (1968) states the accuracy specifications 
recommended for hydrographic suiveys. The depth measure- 
ment specifications are listed in Table I. 
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Table I - Depth Measurement Specifications 
Recommended by the International Hydrographic Bureau 



0-20 meters (0-11 fathoms) 
20-100 meters (11-55 fathoms) 
Deeper than 100 meters 



Depth 



Allowable error 
0.3 meters (1 foot) 

1.0 meters (0.5 fathoms) 
1 % of depth 



The Hydrographic Manual of the National Ocean Survey 
(Umbach, 1976) adds that accuracies attained for all hydro- 
graphic surveys conducted b;^r the National Ocean Survey shall 
equal or exceed the specifications recommended by the Inter- 
national Hydrographic Bureau. These standards do not 
necessarily apply directly to the accuracy requirements 
for a mathematical model of the bottom, but they are good 
reference figiires. 

Solution of the dense data problem for nautical charting 
was the primary motivation for the investigation, but there 
are other uses for models approximating ocean bottom topo- 
graphy. Many coastal processes are closely related to 
bottom topography. These include wave height, wave refrac- 
tion, energy dissipation, wave runup, storm surge and 
beach erosion. Design of offshore structures requires 
input of bottom characteristics. Subsurface, as well as 
surface navigation, could be aided by an ocean bottom model 
stored in an onboard computer. The acctiracy requirements 
and model scales for these applications would be different 
but the modeling methods could be the same. 
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C. SCOPE OF V/ORK 



There are several ways to represent s\irfaces by mathe- 
matical expressions. Those that seemed most applicable 
to the problem are discussed in Section II. Three of the 
models were chosen for experimental analysis. Portions of 
four hydrographic siurveys conducted by the National Ocean 
Survey were used as experimental data sets for this analysis. 
These data sets represent a variation from extreme bottom 
relief to a very flat bottom. The models developed for 
these areas v/ere analyzed quantitatively by comparing 
observed survey depths and computed model depths at the same 
location. Qualitative comparisons of depth contours from 
the two sources were also made. For each type of model, 
the input parameters were varied to investigate minimum 
requirements for a good representation, 

Deteimining the exact location of the shoreline and 
other boundaries is an important part of any survey, but 
including this in the models is beyond the scope of this 
investigation. All the areas used for experimentation 
v/ere restricted so that they do not include shoreline. 
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II. SURPACE MODELING METHODS 



Analytical expressions have been used previously to 
approximate topographic surfaces. Some techniques used in 
map analysis are also applicable to the problem and there 
are some appealing methods which have been used for other 
surface approximations but not for terrain models. None 
of these methods have been used to represent hydrographic 
surveys. Ocean bottom topography is often similar to land 
topography but the research on terrain models has generally 
been for small scale large area maps. The large scale 
hydrographic surveys which must represent detail on the 
order of tenths of fathoms or feet are quite different than 
those large area maps, so modeling techniques which are good 
for small scale terrain models may not be appropriate for 
hydrographic survey modeling. Some important properties of 
the methods v/hich must be considered aside from accuracy are: 

• ease of computation - Must a large system of equations 

be solved to develop the model? 

• dependence of horizontal scale - Hydrographic surveys 
and marine charts of different scales often overlap or 
are adjacent. Eor this reason, it is not good if the 
accuracy of a modeling method varies with horizontal 
distance scale. 

• global versus local models - A global model represents 

a large area with a single expression. A local modeling 
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method represents many adjacent small areas with, 
many corresponding expressions. Generally, there 
is more computation involved in global methods, whereas, 
local modeling reqxiires more data searching to find 
the appropriate local parameters. Global models 
which attain significant data storage savings are of 
particular interest in this study, 

• interpolation versus approximation - Interpolation 
methods generate a surface which fits some data points 
exactly and is used to interpolate between those points 
for surface values at other positions. Approximation 
methods generate a surface v/hich approximates all the 
data but may not fit any data points exactly, A "best 
fit" by some criteria such as least squares is usually 
used. Approximation methods may not represent the 
least depth in an area accurately or they may move the 
position of peaks and deeps significantly. It is impera- 
tive that the model can be controlled to represent criti- 
cal data points exactly. Interpolation methods are thus 
more appropriate for this application. The data points 
which are selected for interpolation will be called 
model points in this presentation. Quite often they 
are significant data points such as a least depth or 
an area of slope change. 

The follov/ing sections discuss methods and previous 
research which are applicable to the problem. 
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A. POLYNOMIALS 



Czegledy (1977), Hardy (1971), Krumbein (1966), and 
V/hitten (1970), discuss the use of polynomials for surface 
representation, A polynomial mapping equation of two inde- 
pendent variables with a specified degree can be produced 
which fits a few data points exactly or approximates all 
the data in a least squares sense. In either case, the sys- 
tem of equations v/hich must be solved becomes ill-conditioned 
as the degree of the polynomial increases. This can be 
alleviated by using orthogonal polynomials. In the method 
of orthogonal polynomials, a collocated series of inde- 
pendent surfaces, linear, quadratic, cubic, etc., is generated. 
The summation of these surfaces is the mapping equation v/hich 
defines the model. Increasing detail is gained by solving 
for and adding the surface of next higher order. This method 
has proven useful for trend analysis of maps. However, it 
has been rejected by some investigators for applications 
requiring more accuracy. The reason as stated by Hardy (1971) 
is that the "ordinary collocated polynomial series is 
unmanageable in representing the sometimes rapid and sharp 
variations of real topographic surfaces." Requiring a high 
degree polynomial to fit closely spaced irregular surface 
pcdits in one area causes significant invalid variations in 
other areas. To avoid these problems, low degree poly- 
nomials have been used in a local approximation mode with 
success, but this does not produce a global surface model. 
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B. DOUBLE EOURIER SERIES 

The double Fourier series model is discussed by James 
(1966) and Krumbein (1966). It is produced by a series of 
independent harmonic surfaces having wave forms of diminishing 
wave length as the order of the sxirface increases. This 
technique has proven valuable for trend analysis particularly 
when the surface features show oscillating patterns. Unfor- 
tunately, the models require high order surfaces to repre- 
sent sharp terrain features. Such surfaces produce oscillations 
v/ith large variations between data points and have many of 
the same drawbacks as the collocated polynomial series. 

i 

C. FINITE ELET'-IENTS 

Gold, Charters and Ramsden (1976) discuss a method of 
surface representation in vmich a system of triangles with 
data points at the vertices is imposed on a surface. An 
interpolating function is used to estimate the surface in 
each triangular element. The interpolant is developed so 
that the stirface passes through the vertices and makes a 
smooth transition from one triangle to the next. 

Peucker, Fowler, little and Mark (1977) have developed 
a similar system of surface representation by Triangulated 
Irregular Networks (TIN). Rather than a smooth interpolant, 
the TIN system uses the planes defined by the three fimction 
values at the vertices of each triangle to represent the 
surface. Considerable work has been done on automated 
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techniques for selecting appropriate points to be used 
for vertices and on development of data structures for 
storage of the vertices, neighboring points, and neigh- 
boring triangles. The TUT system was developed specifically 
for digital representation of topographic surfaces. 

Finite element systems such as these are local methods. 
Detail can be easily incorporated into the model by adding 
points where required v/ithout affecting the model elsewhere. 
Yexj little computation is required but searching the data 
structure to find the appropriate element is necessary. 

Such systems are generally independent of scale unless a 
scale dependent interpolant is used. A single expression 
which represents the surface is not generated by these 
methods. 



D. SHEPAED’S FORIOTIA 

Shepard's method as described by Ppeppelmeier (1975), 
Barnhill (1977) and Franke (1979), has been widely used to 
interpolate random data but has never been used for topo- 
graphic surface representation. The model is produced 
by taking a weighted average of the model points to inter- 
polate the surface value at other points. 

Shepard's formula is expressed by 



f = 




if dj^ = 0 for any i 



if dj^ 0 for all i 



( 1 ) 
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where the axe the depths at the model points; d^^ is the 
distance from the ith model point to the point of computa- 
tion; and the weight assigned to each model point, w^, is 

a function of Two such weighting functions used in 

^i 

this project were simply the inverse distance (1/d^) and 
the inverse distance squared (l/d^ ) . 

In this method, all model points contribute to the 
value of f, but the effect of any model point on the inter- 
polant decreases as the distance from that point increases. 
Another appealing feature of this method is that the value 
of f will always be between the minimum and maximum values 
of the model points. 

Pranke and Little's modification to Shepard's method 
restricts the weighted summation to only those model points 
v/ithin a radius R of the computation point. With this modifi- 
cation, the weighting function approaches zero as the dis- 
tance approaches R and remains zero at distances greater than 
R, The modified Shepard's formula is expressed by 



f = 



or 



f = 





(R-di)^ 

T" 



l^d. 



^i 



if d^ 0 for all i 

( 2 ) 

if dj^ = 0 for any i 
if dj^ 0 for all i 



if dj^ = 0 for any i 



( 3 ) 
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where 



R-d 



(R-di)^ = 



^ for < R 



for d^ > R 



(4) 



The weighting fxuictions l/d^ and produce surfaces 

-m: 

with cusps at the model points. The weighting functions 




and 



~~T~T 



produce siirfaces with flat spots at those 



points. For higher order functions of l/d^ these flat spots 

increase in size and the slopes between them become steeper. 
These properties are shown in two dimensions in Figure 1. 
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Figure 1 - Shepard’s Formula with Various 
Weighting Functions 
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These formulas do not require solution of systems of equa- 
tions and are easily modified by simply adding significant 
data points without recomputing any coefficients. They are 
independent of scaling, global in natxire, and the computa- 
tion is very simple. 

E. DUCHON'S IfETHOD 

The method of Duchon (1976) which was developed as thin 
plate surface theory is described by Meinguet (1979) and 
Harder and Desmarais (1972). It has never been used for 
topographic surfaces but has been used for other surface 
analyses. To develop this model, individual surfaces called 
basis or kernel functions, v/hich are centered at the model 
points, are summed to yield a global surface. There is a 
coefficient associated v/ith each kernel function which 
determines the magnitude of the effect of that kernel ftmc- 
tion on the total stirface. 

The expression for the model is 



i _ ^ C. J 

where n is the number of basis functions and model points 
used. The last three terms represent a plane which is also 
added into the model. The n+3 coefficients C^, i=l, ..., n, 
^2 and are determined by solving the following system 
of n+3 equations. 



f = 




+ A^ + A^ + A^ Y ( 5 ) 
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^1 = 

• 


n 


C, [i(Xi, Y^, X,, Y,*^ ^ . 


+ A^Y^ 


• 

^n = 


1 — ( 

II 

•H 


^i ^n» ^i» ^i^ " Ai . 


A^X^ . A3Y^ 


0 = 


n 

i=l 


"i 


(6) 


0 = 


II hti 
H 


CiXi 




0 = 


n 

i=l 


CiYi 




where f^, f2> 


— 


, f are the surface values at 
^ n 


the model 


points (^lYi) 


CM 







Dudion used two "basis functions 



F (X, Y, X,, I,) , ^,3 

and 

P (X, Y, X., Y.) = d^2 log 

where o o i 

d^ =((X-X^)2 + (Y - Y^r) 



(7) 

( 8 ) 



d^ is the horizontal distance from each model point to the 
point of computation. 

Duchon's method using the above basis functions is 
independent of scale. During experimentation, a third basis 
function 

F (X, Y, Xj_, Y^) = d^ log d. (9) 

was also used. The models using this latter basis function 
are dependent on scale. 
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P. HARDY'S MULTI QUADRIC ANALYSIS 

Multiquadric analysis, as discussed by Hardy (1971, 
1972a, 1972b, 1975 and 1977) resoilted from a search for a 
satisfactory and efficient method to represent topography 
by an analytical model. As suggested by its name, the 
method consists of summing many quadric surfaces (cones, 
hyperboloids, paraboloids, etc.), each associated with a 
model point, to obtain a global surface. Superficially, 
this method is similar to Duchon's method except that the 
kernel functions are quadric surfaces and the additioanl 
three terms are not used. The expression for this model 
is 

f = ^ c. [q (X, Y, X^, Y^^ (10) 

v/here f is the suxface value at the point (X, Y); Q is the 

quadric surface or kernel function; (X^, Y^), i=l, , n 

are the model points at which the kernels are centered; 
and C^, i=l, , n are coefficients assigned to each sur- 

face. 

The following system of equations is used to solve for 
the unknown coefficients. 
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The resulting surface will fit the data exactly at the model 

points (X^, Y^, fj_), i=l, , n. 

Hardy (1971) found that the quadrics which yield the 
best reults are hyperboloids, cones and inverse hyperboloids. 
A hyperboloid is represented by 

Q (X, Y, X^, Y^) = ((X-Xj^)2 + (Y-Y.)2 + 5 ^)^. (12) 

Cones are special cases of hyperboloids where (5 i^ zero: 

Q(X, Y, X^, Y.) = ((X-X^)2 + (Y-Y.)2)'^' = d^. (13) 

dj^ is the distance from the point of computation (X, Y) to 
the center point of the each quadric surface (X^, ^i^* 

Inverse hyperboloid kernels are expressed by; 

Q (X, Y, X^, Y^) = ((X-X^)2 + (Y-Y^)2 + §2)-!^ 

The magnitude of the coefficient determines the steepness 
of the cone or hyperboloid. The sign of determines 
whether the surface is oriented upv/ard or downv/ard. The 
magnitude of 5 determines how flat the hyperboloid is at 
its center. These properties are shown in Figure 2. 
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Figure 2 - Hsrperboloid Kernels 
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For an inverse iiyperlioloid 0 determines tlie peakedness of 
tlie surface at its center point (Figure 5). simply 

represents a multiplicative constant v/hich scales th.e size 
of the surface and specifies its orientation upward or 
dov/nward . 



Q 




Figure 3 - Inverse Hyperboloid Kernels 

For all and d the inverse hyperboloid approaches zero 
as the distance increases. There is an inflection point 
at d = 

The way several quadric surfaces sum to form the global 
surface can be seen in two dimensions in Figure 4. 
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Figure 4 - Quadric Summations (hyperboloid kernels) 



Multiquadric modeling is independent of horizontal 
scaling so long as 5 is linearly related to the horizontal 
distance scales. 
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III. HYDHOGRAPHIO SURVEYS 



The data used in this project were from hydrographic 
sxirveys conducted in the late 1970's hy the National Ocean 
Survey. The methods and procedures used were ty-pical auto- 
mated survey procedures as docvunented hy Umhach (1976) and 
V/allace (1971). A brief description of these procedures 
followed hy more specific information on each data set 
follows. 

A. GENERAL- DESCRIPTION 

Safety of navigation is the primary purpose for which 
hydrographic surveys are accomplished. The data is acqiiired 
hy running sounding vessels in parallel or nearly parallel 
tracklines on the ocean surface and taking depth measure- 
ments along these lines at evely spaced intervals. Cross- 
lines are run at large angles to the main system of lines 
as a gross check on the validity of the data. V/hen indica- 
tions of critical bottom detail are found, development lines 
are run at closer spacing to determine the least depth and 
verify the nature of the feature. The depths are plotted 
on a survey sheet, a sample of which is shown in Figure 5. 

There are many properties of a survey which affect its 
usefiilness. Those most important to this project are s^lrvey 
scale, horizontal positioning accuracy and depth accuracy. 
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Figure 5 - Portion of a Hydrographic Survey Sheet 
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1, Survey Scale 

The survey scale is the ratio of distance on the 
survey sheet to the corresponding distance on the earth. 

The scale chosen for a survey depends on "the area to be 
covered and the amount of detail necessary to depict ade- 
quately the bottom topography and portray the least depths 
over critical features." (Umbach, 1976) The survey scale 
is usually at least twice as large as the scale of any 
chart published for the area, large scale surveys cover 
less area than small scale surveys but greater detail can 
be represented. For this reason, large scale surveys are 
conducted in harbors, anchorages, restricted navigable 
waterways, and areas where dangers to navigation are 
numerous. Areas with considerable detail are 'the most 
difficult to adequately represent by a mathematical expres- 
sion. Three of the four data sets used in this project 
were from large scale surveys. 

2. Horizontal Position Accuracy 

Umbach (1976) specifies that plotted positions, 
"whether observed by visual or electronic methods, combined 
with plotting error shall seldom exceed 1.5 mm (0.05 in.) 
at the scale of the survey." On a 1:5000 scale survey, 
the position of each sounding should thus be represented 
to within 7.5 meters of its actual position on the earth. 
This is important in evaluating a mathematical model. One 
of the data sets had some very steep slopes, where an error 
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of a few meters in positioning would produce a depth 
variation of several fathoms, For areas such as these, a 
much greater depth discrepancy between the model and the 
survey data should be tolerable, 

3, Depth Accuracy 

As seen in Figure 6, there are many components 
that make up the depths represented in a hydrographic 
survey. In addition to the depth recorded by the sounding 
instrument, there are corrections for velocity of sound in 
the water coliunn, the stage of the tide, and the dynamic 
vessel draft. Sometimes surveys have slight inconsistencies 
where data from two different vessels or two different days 
are adjacent or intermixed. These might be due to changes 
in the water column structure that affect the velocity of 
sound, an error in determining offshore tide corrections 
from tide gages near the shore, unrecorded changes in vessel 
speed affecting the dynamic vessel draft, a slight systematic 
error in vessel positioning, etc. Even more critical is 
the effect of waves on the sounding vessel. Small vessels 

/ 

change vertical position rapidly as waves pass while the 
instruments record the depth of the water colirnn below the 
vessel. This depth is too great if the vessel is on a 
wave crest and too small if the vessel is in a trough. The 
angular orientation of the vessel is also affected by waves. 
If the vessel rolls to an angle greater than the sounding 
beam width, the depth recorded may not be xmder the vessel 
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Actual depth 




Pigtire 6 - Components of a Depth Measurement 
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but off to the side. In order to correct for this, the 
echograms were manually scanned, wave action was visually 
meaned out of the record, and depths that were automatically 
acquired with an error greater than the recording interval 
were rerecorded. On days of moderate to heavy wave action, 
this procedure leads to an inordinate amount of manual 
work introduced into an otherwise automated system. Table II 
indicates the depth recording and correction intervals used 
by NOS, Note that in many cases soundings from 0-20 fathoms 
need only be recorded to the nearest v/hole foot or nearest 
half of a fathom. 

Although depth measurement errors can exist on all 
surveys, they are more apparent in areas of flat regular 
bottom. If two adjacent soundings each have nearly a foot 
of error of opposite sign, this will appear as a sharp dis- 
continuity on a flat bottom whereas it will hardly be 
noticed on a steep slope. For this reason, models for 
areas of flat regular terrain when compared with the survey 
data may show some relatively large differences due to the 
data acquisition procedures. 

B. DATA SETS 

Pour data sets were used for analysis in this project. 
They were specifically selected for the variety of bottom 
topography which they represent. 
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TABLE II - Depth Recording and Correction Intervals 
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1. Monterey Bay« California 



This data set was taken from survey registry number 
H-9808. It was conducted in 1979 by NOAA Ship DAVIDSON and 
Naval Postgraduate School personnel and equipment. It 

covers the southernmost part of Monterey Bay including 

« 

Monterey Harbor. The survey was conducted at a scale of 
1:5000. Only one vessel was used on the portion of the 
survey chosen for analysis. The sounding units are fathoms 
and depths range from 0 to 16 fathoms. The bottom has a 
large amoimt of detail. It slopes moderately downward from 
the shore and consists generally of mud and sand. In the 
middle there is an area thick v/ith kelp which is attached 
to a rocky irregular bottom. There are a few rocky areas 
in the deeper part as well. Pigure 7 shows the bottom con- 
tours in one fathom increments. The scale of the plot has 
been reduced for presentation herein. 

2. Morro Bay, California 

This data set was taken from survey registry number 
H-9737. It was conducted in 1978 by the NOAA Ship PAIRV/EATHER. 
It covers a small part of Morro Bay and some navigable water- 
ways open to the bay. The survey was conducted at a scale 
of 1:5000. The sounding units are feet and depths range from 
16 to 82 feet in the portion used for analysis. Pigure 8 
(reduced scale) sliov/s the bottom contours in three foot incre- 
ments. There is one major feature near the center and con- 
siderable irregularity in the northeast comer of the area. 
Otherwise, the bottom slopes gently offshore. 
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Pigure 7 - Monterey Bay Data Set Contours (fathoms) 
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Figure 8 - Morro Bay Bata Set Conto\irs (feet) 
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3. Auke Bay, Alaska 

This data set came from a thesis project by Seidel 
(1979), a student at the Naval Postgraduate School, which 
investigated the affects of using multiple sounding beam 
widths for hydrographic surveys. The procedures were some- 
what non-standard since sounding lines were run much closer 
than normal in an attempt to gain 100% bottom coverage. 
Specifications for 1:5000 scale surveys were used but due 
to the dense sounding spacing, it was plotted at a scale 
of 1:2500. The data was incorporated into survey registry 
number H-9818. It v;as conducted in 1979 by Seidel and the 
NOAA Ship RANIER. It covers a small portion of Auke Bay in 
southeast Alaska. The sounding units are fathoms and depths 
range from 0 to 24 fathoms. The bottom is mostly mud and 
rock and shows a tremendous amount of variation due to 
glacial action. Very steep slopes are encountered in the 
area. At one point, the depth changes from 7 to 22 fathoms 
in a horizontal position change of only 30 meters. Figure 9 
(reduced scale) shows the bottom contotirs of the central 
part of the data set in one fathom increments. 

4. &ulf Coast 

The fourth data set was taken from sxirvey registry 
number H-9785. It was conducted, in 1978 by the NOAA Ship 
MT MITCHEII at a scale of 1:20000 and covers an area in the 
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Gulf of Mexico off the coast of Louisiana. The sounding 
units are feet and depths range from 29 to 37 feet in the 
portion used for analysis. Figure 10 (reduced scale) shows 
the "bottom contours in one foot increments. The bottom is 
generally flat with a very gentle slope. It consists 
mostly of mud and shell fragments. Some of the irregulari- 
ties seen in the bottom contours are in areas v/here the work 
of two vessels overlapped because of crosslines or junctions. 
The flat bottom and small contour increment make these 
irregularities stand out. The survey party reported that 
v/ave action was also a considerable problem during the con- 
duct of this survey. 
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Figure 9 - Auke Bay Data Set Contours 




Figure 10 - Gulf Coast Data Set Contours (feet) 
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IV, RESEARCH PROCEDURES 



A. COICPUTER SYSTE^l 

All the computer work for this project v/as done on 
the IBM 360/67 system at the Naval Postgraduate School's 
W. R, Chxirch Computer Center. All programs were written 
in FORTRAN IV. The VERSATEC-07 electrostatic plotter was 
used for all the data and contour plotting. Both IMSl 
(International Mathematical and Statistical library) 
routines and other library routines v;ere used in the pro- 
grams, 

B. DATA SET PREPARATION 

1. Original Data Condition 

All four data sets were supplied by the National 
Ocean Survey on non-labeled unblocked magnetic tapes in the 
NOS standard record format. Positions of all soundings were 
given in terms of latitude and longitude. Corrected soundings 
were supplied to the nearest tenth of feet or fathoms. Each 
data point had a record sequence number assigned. The NOS 
format also included original observed data and all correc- 
tions to it as v;ell as descriptive cartographic codes and 
other information. 
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2, Program TAPCNV - Tape Conversion 

Only corrected position, corrected depth and 
record number identification were required for this pro- 
ject. The program TAPCNV was written to read this data 
from the non-labeled NOS tapes. The geodetic positions 
were converted to an X-Y plane coordinate system based on 
the Modified Transverse Mercatot (MTM) projection (Wallace, 
1971). Double precision computations were used for this 
conversion. The MTM projection gives the positions in 
terms of meters of northing and meters of easting from 
a local origin. This X, Y position was then converted to 
plotter coordinates in terms of inches from the plotter 
origin. The record mmber, depth, geodetic position, 

MTI^ coordinates and plotter coordinates were blocked and 
recorded on disk and on an NPS tape with standard system 
labels, 

3. Program DATPIT - Data Plotting 

This program was written to display the discrete 
point data on a plotted sheet. latitude and longitude 
grid intersections at specified intervals were converted 
to plotter coordinates and straight lines v/ere drawn con- 
necting these points to provide the geodetic position 
reference system. Two sheets were plotted with this 
reference grid. On one sheet depths were plotted to the 
nearest tenth, (NOS plots tenths only in shallow water 
v/hen the depth units are fathoma) Record numbers were 
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plotted on the second sheet. Overlaying the two sheets 
facilitated reference to any particular data point. DATPIT 
v;as used to plot the entire siirveys as a first step. V/hen 
portions of the surveys were selected for analysis DATPIT 
was used again to select and plot only the data within 
the specified area. 

4. Program CONPAT - Data Contouring 

Part of the data anal.ysis consisted of comparing 
contours of the original data with contours from the model. 
Initially, contours of the survey data were hand drawn - a 
procedure that is somewhat subjective. In order to remove 
as much subjectivity as possible from the analysis hand 
contouring was replaced by machine contouring. The library 
routine CONISD, for contouring irregularly spaced data, was 
used in the program CONPAT. This routine first generates 
triangles with data points at the vertices. By linearly 
interpolating along the triangle sides for the contour 
values, points on each contour are found and connected to 
generate the contour lines. The contotirs generated in this 
manner were not smooth as would be desirable, but the data 
points were dense enough so that this was not a problem. 

C. MODEL DEVELOPMENT AND ANALYSIS 

Program MODEL was written to do the model development, 
the quantitative analysis, and to aid in the qualitative 
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analysis. To change from one modeling technique to another 
the only modification necessary was the replacement of one 
module. That module contained the routine to develop the 
model by a given method and to compute the depth at any 
point. All the modeling techniques required the selection 
and use of model points from the survey data. The model 
points were specified by record numbers on punched card 
input and the survey data was read from disk and stored in 
memory. The model points were stored in arrays for model 
development , 

1. Coefficient Computation - Subroutine LEQ2S 

The methods of Hardy and Buchon require solution 
of S3rmmetric systems of linear equations to determine the 
model coefficients. The double precision version of the 
IMSL routine 1EQ2S was used for this purpose. This routine 
uses syimnetric decomposition with iterative improvement 
to solve the systems. Systems of up to 226 equations in 
226 unknowns were solved during the course of this project. 
The model coefficients and respective model points v/ere 
output for analysis. 

2, Quantitative Analysis - Subroutine STAT 

Subroutine STAT was developed to provide a quanti- 
tative analysis of each model. Each survey depth was com- 
pared with the depth computed from the model at the same 
X, Y position. The root mean square difference, maximum 
positive difference and maximum negative difference were 
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tabulated for each run. A positive difference signified 
that the model depth was deeper than the survey depth; a 
negative difference signified that the survey depth v;as 
deeper. The root mean square difference is given by the 
expression 



v/here IlD is the model depth, SD is the survey depth, and 
N is the number of points used for the comparison. Those 
points with a difference greater than 1,5 times the RMS 
difference were listed for manual inspection and analysis. 

3. Qualitative Analysis - Subroutines SETCON and COITTTJR 
The qualitative analysis was accomplished by compari- 
son of model contours with those from the original data. 

In order to produce the model contours subroutine SSTCON 
developed a quarter inch grid over the modeled area at the 
scale of the survey. The model depth was computed at each 
grid intersection. These depths and positions were passed 
to the library routine CONTUR for contouring. CONTUR is 
similar to COniSD except that it was written specifically 
for gridded data and runs considerably faster than CONISD. 



RMS difference 




(15) 
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Y. RESEARCH RESULTS 



A. SELECTION OE IffiTHODS FOR EXPERIMENTATION 

Of the modeling techniques discussed in Section II, 
Duchon's method, Shepard's formula and Hardy's multiquadric 
analysis were selected for testing. Neither Euchon's method 
nor Shepard's foimiula had previously been used for terrain 
surfaces. Multiquadric analysis had been used with good 
results for topographic data but not at the scales and 
accuracy requirements necessary for hydrographic survey 
representation. 

The methods of polynomials and double Fourier series 
were not tested. They had proved to be useful for some 
applications such as trend analysis and representation of 
repeated features. The fact that forcing polynomials or 
double Fourier series to fit irregular data in small areas 
produces unwanted irregularity in other areas would seem 
to preclude these methods from producing good results in 
this application. Some methods reduce this effect by using 
local expressions which fit only small areas at a time, but 
this defeats the purpose of generating a global model to 
represent large areas of the data. 

Finite element methods were not tested. They are strictly 
local methods which, in addition to storage of model points, 
reqxiire storage of pointers to the neighboring points and 
neighboring triangles of each model point. 
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B. I4ETH0D COl-rP ARISON PROCEDURES 



The Monterey Bay and Morro Bay data sets were used for 
comparison of the methods. The Monterey Bay data set was 
used in the qualitative manner only. The Morro Bay data 
set was used for Both qualitative and quantitative compari- 
son. The procedures described below v/ere used for all methods 
in order to make controlled comparisons. 

As a first step, 42 data points on a 6x7 grid v/ere chosen 
from the Monterey Bay data at regular spacing v/ithout regard 
to bottom detail. After the models were generated with these 
points, an additional 18 model points were selected in areas 
where more detail v/as required. Models were then generated 
using the 60 points. The third step was to choose 30 more 
points aroimd the outside of the original area at the same 
spacing as the original 42 points extending the grid to 8x9. 
Models v/ere generated with the 90 points to determine the 
affect of extending the model area. 

Thirty-seven model points were chosen at regular spacing 
from the Morro Bay data set in the first step with that 
data. Thirty and thirty-one additional points were selected 
in the second and third steps for totals of 67 and 98 model 
points. The additional points were all within the original 
area in places v/here additional detadl or accuracy was needed. 
The effect of increasing the model point density was examined 
in this way. The statistical results, as well as contours 
of the Morro Bay tests, were compared. 
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C. RESULTS OF DUCHON'S METHOD 



Equation 5 shows that Duchon's model is produced by the 
summation of a plane and a series of basis fionctions centered 

at the model points, Duchon's method with the basis functions 

2 

d and d logd and a similar method with basis fimction d logd 
were tested, 

1, General Findings 

For all three basis fxmctions the 42 Monterey Bay 

points produced similar models which showed the general 

trend of the bottom but Tery little detail. Using the 60 

3 2 

model points, the basis fimctions d*^ and d logd gave more 
detail and a fair representation of the bottom trends. 

See Figure 11, 

The basis fiuiction d logd gave a very poor repre- 
sentation of the bottom when the 60 points v/ere used. In 
an effort to resolve detail in some areas, several points 
v/ere chosen very near each other. This caused steep slopes 
to be generated which extended into areas where there were 
no model points and created invalid peaks and deeps. 

The quantitative results of Duchon's method using 
the Morro Bay data set are given in Table III, Note that 
the model using the basis fvmction d logd became better in 
both RI4S difference and maximum difference as the number 

of points was increased. The results didn't change much 

2 

using the basis fiuiction d logd and they became worse for 
the basis function d , The contours reflected these 
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3 2 

statistics. In the cases of the d and d logd basis fxmc- 
tions, the detail was increased in the areas where points 
were added but the statistical results did not improve due 
to greater error in other areas. As shown in Pigure 12, 
the model using 98 points with the basis function d logd 
compares very favorably with the original data contours in 
Pigure 8, 

TABLE III - Results of Duchon’s Method - Morro Bay 



Basis 
Pune ti on 


Nr of 
model 
points 


PUS 

difference 


Maximum 

positive 

difference 


Maximum 

negative 

difference 


Nr of 

data 

points 


d3 


37 


1.20 


3.89 


-6.15 


936 


d^logd 


37 


1.14 


3.73 


-6.29 


936 


d logd 


37 


1.15 


4.32 


-6.51 


936 


• d5 


67 


1.77 


6.74 


-7.99 


936 


d^logd 


67 


1.11 


4.09 


-4.42 


936 


d logd 


67 


0.77 


3.11 


-3.15 


936 


d^ 


98 


2.13 


7.75 


-17.83 


936 


d^logd 


98 


1.14 


4.19 


-6.15 


936 


d logd 


98 


0.67 


2.23 


-2.48 


936 


2. 


Dependence 


1 on Scale 









The results of the previous section lead to a ques- 
tion concerning the ability of the basis function d logd to 
produce a much better model than other basis functions in 
one case but not in the other. The reason for this turned 
out to be the scale of the data. The first tv/o basis 
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Pigvire 12 - 98 Point Model of Morro Bay 
(basis function d log d) 
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functions produce models which are independent of scale 

whereas, the third produces models which are not. Initially, 

the Morro Bay data was used at the earth’s natural scale in 

meters. At this scale, the system of equations was ill- 

conditioned to such a degree that it couldn’t he solved. 

The data presented in the previous section was acquired 

with the horizontal position data scaled to a distance of 

unity on a diagonal from one comer of the area to the 

opposite comer. In this case, the results were good. 

The Monterey Bay data was scaled for a diagonal distance 

of 50, The results in this case were poor. Table lY gives 

the results of some tests nm at various scales using the 

2 

Morro Bay data and the basis functions d logd and d logd. 

The set of 98 model points and 956 data comparison points 
were used for these tests. 



58 



TABLE IV - Effects of Scale on Duclion's Method - 

Morro Bay- 



Basis fimction: 


d logd 






Diagonal 

distance 


RMS 

difference 


Maximum 

positive 

difference 


Maximum 

negative 

difference 


0.001 


0.675 


2.115 


-2.365 


0.01 


0.672 


2.121 


-2.379 


0.1 


0.668 


2.131 


-2.406 


1.0 


0.671 


2.216 


-2.479 


10.0 


5.829 


13.566 


-20.285 


100.0 


17.394 


65.291 


-198.561 


1000.0 


0.781 


2.211 


-2.461 


Basis f-unction; 


d^logd 






0.001 


1.138 


4.186 


-6.149 


0.01 


1.137 


4.185 


-6.143 


0.1 


1.139 


4.193 


-6.155 


1.0 


1.138 


4.189 


-6.151 


10.0 


1.138 


4.189 


-6.151 


100.0 


1.138 


4.189 


-6.151 


1000.0 


1.138 


4.190 


-6.153 
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D. RESULTS OP SHEPARD'S METHOD 

As defined in equation 1, Shepard's formula is a 
weighted summation of the model point depths. The weight 
assigned to each model point is a function of the inverse 
distance from the point of computation to the model point. 
The weighting functions used in this analysis were: 



where (R-dj^)_^ is defined in equation 4. 

1. Computation of R 

In the modified Shepard's method, a radius of in- 
fluence, R, is used. Rather than choosing this radius 
arbitrarily, a method was used which related R to the den- 
sity of the model points independent of the scale of the 
data. This also allowed variation of R according to the 
average number of model points which v/ould fall within the 
radius of influence. 



The following expression was used for this purpose: 




1 




R 




(16) 
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where D is the maximum distance between any two model 
points, N is the total number of model points, and NPPR 
is the average number of points which should fall within 
the radius of influence, related to the size 

of the area v/hich contains the data points. Dividing 
this by N gives a measure of the average area which could 
be assigned to each point. Multiplying by NPPR gives the 
area which could be associated with that many data points, 

T a.king the square root of this gives a radius v/hich would 
define that amomt of area centered at the point of compu- 
tation, On the average there should be NPPR model points 
within a distance R from any point of evaluation. The 
tabulated statistical results express the radius in terms 
of NPPR instead of R, 

2, Inverse Distance Weighting Function 

Table V gives the statistical results of the tests 
using the inverse distance weighting functions. The table 
shows that use of the modified Shepard method improved the 
results considerably. In all cases, the best results were 
obtained by including an average of six model points in 
the radius of influence. The table also indicates that no 
statistical improvement was made by increasing from 37 to 
98 model points. 

The contours produced by this method (Pigure 13) 
for both data sets were poor. The basic trend of the bottom 
can hardly be seen. The contours are quite wavy where they 
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TABLE V - Shepard’s Formula with Inverse 
Distance Weighting Function - Morro Bay 



Number Maximum Maximum Number 



of model MS positive negative of data 



points 


NPPR 


difference 


difference 


difference 


points 


37 


6 


2.06 


8.17 


-9.37 


936 


37 


9 


2.28 


7.33 


-9.68 


936 


37 


All* 


10.88 


24.43 


-25.96 


936 


67 


4 


2.45 


8.50 


-7.97 


936 


67 


6 


2.17 


5.69 


-6.94 


936 


67 


9 


2.31 


6.86 


-7.61 


936 


67 


25 


3.53 


9.61 


-12.68 


936 


67 


All* 


11.41 


22.77 


-28.58 


936 


98 


4 


2.41 


8.11 


-8.95 


936 


98 


6 


2.17 


7.37 


-6.73 


936 


98 


9 


2.22 


8.64 


-7.51 


936 


98 


25 


3.46 


11.38 


-12.51 


936 


98 


56 


5.03 


13.33 


-16.50 


936 


98 


All* 


12.50 


21.33 


-32.22 


936 



* Shepard’s formula - All model points contributed to 

the weighted average. 
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PigTire 13 - 98 Point Sheoard's Formula Model of Morro Bay 

(NPPR = 6) 
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should he straight. In some cases, peaks or deeps are pro- 
duced at the positions of model points which aren't found 
in the original data, 

3. Inverse Distance Squared Vfeighting Function 

Table VI gives the statistical results of similar 
tests using the inverse distance squared weighting function. 
There is considerably less variability as NPPR is changed 
using this weighting function. The results are better for 
large NPPR and for the unmodified version, but the best 
results at smaller NPPR did not improve. 

S. RESULTS OP HARDY'S MUITIQUADRIC ANALYSIS 

As indicated in equation 10, Hardy's multiquadric model 
is generated by summing quadric kernel surfaces, each of 
which are centered at model points. Hjrperboloids, cones 
and inverse hj^perboloids were the kernel surfaces tested, 

1. Determination of (5 

Both hyperboloids and inverse hyperboloids require 
the parameter (5 (Section II. P). Variation of (5 makes 
considerable difference in the results. The effect of 
any value of (5 on the shape of the quadric surfaces with 
respect to the entire model is related to the scale of ^ 
the model. Hardy (1977) has indicated that the optimum 
value of 8 in his investigations was also related to the 
distance between model points. The following expression 
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TABLE VI - Shepard’s Formula with Inverse 
Distance Squared V/eighting Function - 
Morro Bay 



Number 
of model 
points 


NPPR 


RMS 

difference 


Maximum 

positive 

difference 


Maximum 

negative 

difference 


Number 
of data 
points 


37 


4 


2.83 


9.51 


-8.90 


936 


37 


9 


2.37 


8.75 


-9.09 


936 


37 


25 


2.36 


8.38 


-9.57 


936 


37 


All* 


5.00 


14.62 


-17.06 


936 


67 


4 


2.81 


9.14 


-8.87 


936 


67 


9 


2.31 


6.01 


-6.88 


936 


67 


25 


2.34 


6.71 


-9.00 


936 


67 


All* 


5.54 


15.42 


-19.90 


936 


98 


4 


2.58 


9.33 


-9,46 


936 


98 


6 


2.33 


7.53 


-7.84 


936 


98 


9 


2.18 


7.75 


-7.05 


936 


98 


25 


2.26 


9.56 


-8.40 


936 


98 


56 


2.72 


10.57 


-11.66 


936 


98 


All* 


6.58 


15.50 


-22.48 


936 



* Shepard’s formula - All model points contributed to 

the weighted average. 
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was used to relate 8 to the average density of the model 



points : 



5 




X 0.1 X NPPR 



(17) 



V/ith this expression, the effect of 8 on models of 

different scales will be similar as long as ITPPR is the 

> 

same. The tables in the following sections are expressed 
in terms of NPPR instead of the absolute value of 8 • 

A cone is a special case of hyperboloid where 8 is 
zero. The tables for hyperboloid kernels include cones by 
listing NPPR as zero. A zero value of 8 is not valid 
for inverse h 3 rperboloids since the peak of an inverse 
h 3 rperboloid increases to infinity as 8 approaches zero. 



2. Inverse Hyperboloid Kernels 

For both data sets when 8 was small (NPPR=5), the 
contours showed holes at the model points v/hich were not 
indicated in the original data. The representation of the 
actual surface was very poor. Increasing NPPR to 10, 15 
and 20 gave somewhat better results and the bottom trends 
v/ere evident but the representation was still not good. 
V/ith NPPR greater than 20, very steep slopes were created 
in large areas where no model points v;ere chosen. 

The statistical results using the Morro Bay data 
set are given in Table VII. The results became worse as 
more model points were added. The best results were con- 
siderably poorer than the best results from other methods, 
particularly the maximum differences. 
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TABLE VII - Multiquadric Analysis with. 
Inverse Hyperboloid Kernels - Morro Bay 



Number 
of model 
points 


OTPR 


MS 

difference 


Maximum 

positive 

difference 


Maximum 

negative 

difference 


Number 
of data 
points 


37 


5 


4.99 


4.71 


-33.48 


936 


37 


10 


2.05 


4.82 


-21.09 


936 


37 


15 


1.59 


4.56 


-14.73 


936 


37 


20 


1.43 


4.24 


-10.90 


936 


67 


5 


6.11 


5.50 


-32.06 


936 


67 


10 


2.54 


11.68 


-21.15 


936 


67 


15 


2.69 


17.74 


-15.71 


936 


67 


20 


4.09 


20.16 


-15.60 


936 


98 


5 


23.85 


2.35 


-60.04 


936 


98 


10 


3.13 


9.92 


-23.96 


936 


98 


15 


3.98 


21.97 


-20.13 


936 


98 


20 


8.23 


67.37 


-33.72 


936 
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3, Hyperboloid and Conic Kernels 



HyTperbolic and conic kernels were evaluated on 
both data sets and NPPR was varied from zero to 25 for 
each set of model points. With 42 regularly spaced model 
points on* the Monterey Bay data set, not much detail was 
evident but the general bottom trends were well represented 
for all values of NPPR. Increasing to 60 model points 
gave more variation with NPPR. For NPPR set to zero and 
one the results were very good. See Figure 14. The 
detail was improved and the bottom trends v/ere still accu- 
rate. For NPPR set to 10 and 20, the results became pro- 
gressively worse. Very steep slopes were generated which 
created invalid peaks and deeps in areas where no model 
points were chosen. The reason for these slopes is apparent 
when examining the magnitude of the coefficients. For 
l'TPPR=l, the mean coefficient magnitude was 0.33; for !IPPR=20, 
the mean coefficient magnitude was 151.63. V/hen NPPR was 
increased to 25, the system of equations became so ill- 
conditioned that it could not be solved. This is due to 
the increased flatness of the h 3 rperboloids v/hen NPPR becomes 
large. In areas where several model points are very close 
in order to represent sharp irregularities in the bottom, 
the flat hyperboloids centered at those points can't produce 
the detail required. Increasing to 90 model points produced 
similar results. 

The statistical results from the Morro Bay data set 
are given in Table VIII. For small NPPR, the results became 
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Figure 14 - 60 Point Multiquadric Model of Monterey Bay 

(hypertoloid kernels) 



TASIE VIII - Multiquadric Analysis with 
H 3 rperholoid and Conic Kernels - Morro Bay 



Nmber 
of model 
points 


NPPR 


RMS 

difference 


Maxim-urn 

positive 

difference 


Maximum 

negative 

difference 


Number 
of data 
points 


37 


0 


1.18 


8.37 


-6.59 


936 


37 


1 


1.14 


6.00 


-6.76 


936 


37 


10 


1.18 


3.97 


-6.92 


936 


37 


25 


1.24 


3.90 


-6.17 


936 


67 


0 


0.75 


3.61 


-3.09 


936 


67 


1 


0.78 


3.59 


-3.01 


936 


67 


10 


2.10 


15.41 


-6.79 


936 


67 


25 


12.04 


62.50 


-44.52 


936 


98 


0 


0.65 


1.92 


-2.14 


936 


98 


1 


0.68 


2.06 


-2.12 


936 


98 


10 


2.84 


12.72 


-14.11 


936 


98 


20* 


14.44 


117.96 


-43.61 


936 



* system couldn't be solved for NPPR=25 
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continua21y tetter as more points were added for greater 
detail. As the model points hecame more dense, the test 
results were acquired hy using conic kernels (NPPR=0). 

For the original 57 regularly spaced model points, the 
best results were for small but non-zero NPPR. The con- 
tour comparisons reflected the model quality demonstrated 
by the statistics. 

P. SUTMARY 

A graphical comparison of the statistical results of 
the methods is given in Figure 15. Duchon’s method with 
basis functions d'^ and d logd gave only a fair representa- 
tion of the bottom with regularly spaced model points. 
Additional model points did not improve the results so 
this technique was rejected. The basis function d logd, 
was introduced which gave good results (comparable to the 
miiltiquadric method) in one case and poor results in another. 
This was due to a dependence on the horizontal scale of the 
data. Independence of scaling for hydrographic survey 
modeling is veiy important since soirveys are plotted at 
various scales. The method with basis function d logd is 
unacceptable for this reason. 

Shepard’s formula gave best results in modified form 
v/ith about siz model points in each radius of influence. 

The inverse distance v/eighting function was better than 
the square of the inverse distance. The results were 
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RMS DIFFERENCE (FEET) 




Figure 15 - Comparison* of Modeling Methods 
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considerably worse than those of Duchon's or Hardy's methods. 
Improvement could not be gained by increasing the number of 
model points. For these reasons, Shepard's method was 
unacceptable. 

Hardy's multiquadric analysis v/ith inverse hyperboloids 
gave very poor results. For small ^ boles were produced 
at the model points and the depths between model points 
were not accurate. 

Of the methods tested, only multiquadric analysis with 
conic or sharply pointed hjrperboloid kernels gave results 
which iadi cate d that further tests were warranted. Depic- 
tion of detail is improved by adding more model points with- 
out adversely affecting the model in other areas and the 
method is independent of linear scaling. 
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VI. PURTHER TESTING OF MULTI QUADRIC MALYSIS 
mUE CONIC MD HYPERBOLOID KERNELS 



The results from the previous section showed that 
multi quadric analysis with conic or sharply pointed hyper- 
boloid kernels was the only method tested that could meet 
the requirements of this application. Additional experi- 
mentation was done to determine the best procedures for 
selecting the model points and for joining models together 
at the boundaries. Tests were also run to determine how 
accurately the data sets could be represented v/ith addi- 
tional model points while still saving significant storage 
space. 

A. SELECTION OF MODEL POINTS 

The selection of the data points to be used for the 
modeling is a critical process in the development of the 
multiquadric model. Three models for selection of the 
points were tested using the Auke Bay, Alaska data set. 
All point selection was done manually but consideration 
was given to the difficulty in automating the process. 

1. Regular Spacing Selection 

In this method, data points from the survey were 
chosen at nearly even spacing v/ithout regard to depth. 
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bottom features, contour separation or any other factor. 

To avoid biasing, they were selected from a plot of record 
numbers rather than a plot of depths or depth contours. 
Additional points for more detail and accuracy were chosen 
for subsequent runs maintaining even spacing as much as 
possible ^d-thout considering any factor except the hori- 
zontal distribution. The results of this procedure are 
presented in Table IX. The RI-IS differences were improved 
significantly v/hen the number of model points was increased 
from 53 to 110 but the maximum positive differences were 
not improved. Additional densification of the model points 
produced little improvement in either the RI4S differences 
or the maximum differences. 

2. Iterative Selection 

In the iterative selection process, the results 
of one model were used to eliminate some model points and 
select additional ones to produce a better model. After 
developing the model with the first set of points, the 
comparison of survey data points with the model v/as analyzed. 
Additional model points were selected wherever single 
point comparisons showed the largest differences or in areas 
where several points showed relatively large differences 
of the same sign. Model points which had very small asso- 
ciated coefficients were eliminated. A small coefficient 
indicates that the associated basis function has little 
effect on the model since it remains near zero within the 
modeled area. The model was then computed v;ith the new set 
of points. 
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TABLE IX - Selection of Model Points with Even Spacing 
Number Maximum Maximum Number 



of model RMS positive negative of data 



points 


ITPPR 


difference 


difference 


difference 


points 


53 


0 


2.04 


9.51 


-5.64 


1407 


53 


1 


2.00 


9.45 


-5.78 


1407 


53 


5 


1.91 


8.92 


-6.39 


1407 


53 


10 


1.90 


8.61 


-7.03 


1407 


53* 


15 


1.91 


8.68 


-7.66 


1407 


53 


20 


1.93 


8.68 


-8.10 


1407 


no 


0 


1.37 


9.84 


-4.25 


1407 


no 


1 


1.33 


9.95 


-4.71 


1407 


no 


2 


1.30 


10.04 


-5.02 


1407 


no 


5 


1.26 


10.15 


-5.46 


1407 


no 


7 


1.26 


10.19 


-5.60 


1407 


no 


10 


1. 26 


10.28 


-5.76 


1407 


no 


15 


1.29 


10.46 


-5.97 


1407 


144 


1 


1.25 


9.86 


-4.26 


1407 


144 


3 


1.21 


9.95 


-4.66 


1407 


144 


5 


1.13 


10.07 


-5.41 


1407 


144 


7 


1.11 


10.07 


-5.58 


1407 


144 


10 


1.10 


10.06 


-5.76 


1407 


144 


15 


1.11 


10.10 


-6.01 


1407 
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This procedure could be repeated \mtil the desired 
accuracy was attained, the maximum n\unber of model points 
to be used was reached, or the model accuracy no longer 
improved with further iterations. Por this comparison of 
selection methods, the process was repeated until the number 
of model points was approximately the same as the maximum 
number used in the test of the regular spacing selection 
method. 

Table X shows the results of these tests. In all 
tests, the best results were obtained v/hen NPPR=0 (conic 
kernels). Both RMS and maximum differences improved sig- 
nificantly as the selection process was repeated. Two 
iterations yielded approximately the same number of model 
points as the maximum used in the regular spacing selection 
method. 

Points related to features such as peaks, deeps 
or sharp changes in slope v/ere chosen for the initial set 
of model points in these comparisons. Regular spaced points 
for the initial set v/ere used in other tests with the 
iterative method. The results were good for both methods 
of initial selection. After a few iterations relatively 
fev/ points from the initial set remained so the initial 
point selection method made little difference. 

A comparison of model point selection by regular 
spacing and by iteration is shown in Figure 16. 



77 



TABLE X - Selection of Model Points by Iteration 



Number 
of model 
points 


NPPR 


RMS 

difference 


Maximiun 
positive 
diff er^ce 


Maximum 

negative 

difference 


Niomber 
of data 
points 


82 


0 


1.48 


4.33 


-5.52 


1407 


82 


1 


1.56 


3.85 


-5.50 


1407 


82 


5 


2.36 


8.09 


-9.43 


1407 


107 


0 


1.06 


3.04 


-3.29 


1407 


107 


1 


1.15 


2.81 


-3.32 


1407 


107 


3 


1.50 


4.89 


-4.82 


1407 


152 


0 


0.72 


1.91 


-2.58 


1407 


152 


1 


0.73 


2.03 


-2.72 


1407 
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Figure 16 - Comparison of Model Point Selection Methods 
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3. Complete Selection Tjy Topographic Feature 



Vfliile using the iterative selection process, it 
was found that the additional points were selected where 
there were significant changes of slope or where there 
were large areas vathout any model points. This led to 
an attempt to select all the model points in one step 
"based on the following criteria. 

• Select points at peaks, deeps, ridges and where 
slopes change significantly. 

• Select points to avoid leaving any large areas 
without model points as a result of the first 
criterion. 

A test was done "by choosing 145 points to model 
only half of the Auke Bay data set. The RE-IS difference 
was 0.88 and the maximum difference v/as -3.22. These 
results show that one-shot selection is not nearly as 
good as the iterative method and probably not much better 
than the regular spacing method. 

4. Summary 

The iterative selection method gave by far the 
best statistical results. It also required the most com- 
puter time. It v;ould be adaptable to complete automation 
since the method for point selection has little subjectivity 
involved. 
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The regular spacing method would he easier to auto- 
mate hut it doesn’t give good results when detail is 
required. The method of complete selection hy feature 
doesn’t give significantly better results than does regu- 
lar spacing and the method would he difficult to automate. 

B. MODEL JTJMCTIONS 

Agreement between two or more stirveys v/hich have a 
common boundary or cover a common area is a very important 
check on the quality of the siirveys. Similarly, agreement 
between the models representing the surveys must be main- 
tained to avoid ambiguity. The problem is even more acute 
vidien several models are joined together to represent a 
single survey. It would be desirable to represent an entire 
survey v/ith a single model but in many cases this v^ould be 
difficult due to the large systems of equations which would 
have to be solved to generate the model coefficients. 

Hardy (1971) suggested a simple method of junctioning 
where common points on the boundary were used in adjacent 
models. This woiild assure that the models were in agree- 
ment at these points and if 'chosen at close intervals, the 
differences at intermediate points on the boundary would 
be relatively small. That method is not appropriate for 
this application since the data points are not along 
straight lines which could be used as boxmdaries. Common 
points on irregular boundaries could be used but this v^ould 
complicate model boundary definition and storage. 
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Tv /0 other more appropriate methods were investigated 
for this application. Both use overlapping areas in the 
adjacent models rather than a common boundary. 

The first method is analogous to the method of using 
common points on a bo^mdary because the model points within 
the overlapping area are required to be the same for both 
models. A line in the middle of the overlapping area is 
chosen to delineate the areas of model usage. See Figure 17. 
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Figure 17 - Model Junctions by Overlap 



This method could produce small discontinuities at the 
boundary. 

The second method eliminates the discontinuity completely 
but requires more computation. The model points in the over- 
lapping area are not required to be common in both models. 
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In this area, a weighted sum of the values obtained from 
each model is used. The weight, v;, is determined by the 
Hermite Polynomial 

w = 1 - 3s^ + 2s^ 

where s is the relative distance from the point of compu- 
tation to the overlap area boundary. The value of s varies 
from one at the outer model boundary to zero at the inner 
boundary.'' of the overlap area (see Figure 18). 
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Figure 18 - Model Jimction by Hermite Poljaiomial 



The weight assigned to the model A value at the point of 
computation in Figure 18 is determined by using 



s 




(18) 
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The weight assigned to the model B value at the same point 
is determined hy using 



s = % . (19) 

The sum of the resultant weights is always one, A plot of 
the Hermite polynomial is shown in Figure 19. 




Because the first derivative of this function is zero at 
s=0 and s=l, the transition from one model to an adjacent 
one v/ill be smooth and continuous. 

The Aulce Bay survey was divided into two overlapping 
models as pictured in Figures 17 and 18. Each was modeled 
separately using the iterative method of model point selection. 
The two models were then joined by the two methods dis- 
cussed above. The results are presented in Table XI. 

Both junction methods showed improved results over the 
individual models since some of the largest errors were 
located near the outer boundaries of the models. The results 
with the polynomial method were only slightlj'’ better. Con- 
tour comparisons between the junction methods showed little 
difference. The possible discontinuity at the boundary when 
not using the Heimite polynomial was not apparent in the 
contours. 
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TABLE XI - Results of Model Junction 
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These results have shown that transition from one model to 
another can he done smoothly. The method using the Hermite 
polynomial is only slightly better than the method using 
common points in an overlapping area, VJhen joining models 
without common points, e.g, two different surveys, the 
^ Hermite polynomial method will give a smooth non-ambiguous 
transition from one model to another. 

C. MODEL REPIHEI'OTT 

Multiquadric analysis and the iterative method of model 
point selection were used to refine the models of all four 
data sets. Tables ]CII, XIV and X^/ give the statisti- 

cal results of these tests. Figures 20, 21, 22 and 23 shov/ 
the effect of each iteration on the Rf-IS differences. The 
ntimber of model points is expressed as a percentage of the 
number of data points represented. This is a direct indi- 
cation of the storage savings attained by each model. All 
four figures show a similar trend. Initially, the results 
improve rapidly as the percentage of data points used in 
the model increases. The improvement then tends to level 
off and repeated iterations generate less improvement. 

Contour plots of the final models for each data set 
are shown in Figures 24, 25, 26 and 27. Comparison of 
these v/ith Figures 7, 8, 9 and 10 shov/ the agreement with 
contours of the original data. 
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For the Monterey Bay model, an Rl/IS difference of 0.2 
fathoms and a maximum difference of -0.74 fathoms are good 
considering the irregularity of the bottom. These results 
were obtained using 17.2% of the data points as model 
points. 

For the Morro Bay model, an El^S difference of 0.55 
feet and a maximum difference of -1.58 feet are good consifer- 
ing that the depths range from 16 to 82 feet. This repre- 
sentation was made using only 15.7% of the data points. 

The allowable error specification (Section 1.3) is one foot 
in the shallow end of this range and three feet in the deep 
end. 

For the Auhe Bay data set, an Rl-'IS difference of 0.30 
fathoms and a maximum difference of -0.95 fathoms using 
20.6% of the data points are good considering the steep 
slopes in the area. A horizontal positioning error of a 
few meters (within tolerance for the survey scale) could 
create several fathoms difference in many of the recorded 
depths. 

Even though the range of depths in the Gulf Coast data 
set is small, and Bl^IS difference of 0.30 feet and a maximtim 
difference of -1.16 feet using 9.9% of the data points 
could be considered good. The allowable error (Section I.B) 
for these depths is one foot. There are places in the data 
set where crossline and adjacent sotindings from multiple 
vessels disagree by as much as two feet. The model 
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representation tends to smooth out such discrepancies and 
this smoothing appears as relatively large differences in 
the statistics. 

The Gulf Coast model gave the only case where a large 
5 (l'TPPR=50) gave much "better results than smaller values. 

Only nine model points were used in that case, V/hen more 
model points were added a small NPPR was required for good 
re suits. 
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table XII - Monterey Bay Model Results 
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TABLE XIII - Morro Bay Model Results 
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TABLE XIV - Atilce Bay Model Results 
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West side of data set 
East side of data set 

East and west sides joined by Herraite Polynomial Method 



TABLE XV - Gulf Coast Model Results 
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TABLE XV (cont) 
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Figure 20 - Monterey Bay Model Results 
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Figure 21 - Morro Bay Model Results 
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Figure 22 - Auke Bay Model Results 
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Figure 23 - Gulf Coast Model Results 
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Figure 24 - 226 Point Monterey Bay Model Contours 
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Figure 25 - 144 Point Morro Bay Model Contours 
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VII. CONCLUSIONS 



Of tile methods examined, only Hardy's multiquadric 
modeling technique with conic kernels was found to he suit- 
able for modeling hydrographic survey data. Polynomial, 
double Pourier series and finite element methods were 
rejected for reasons found in the literature. Duchon’s 
method, Shepard’s formula, and multiquadric analysis with 
inverse hjT)erboloid kernels were rejected as a result of 
tests presented herein. 

Selection of model points for the multiquadric method 
is best performed by iteration. An initial set of model 
points may be chosen either by regular spacing or by bottom 
feature. Comparisons of the initial model with the survey 
data are used to select additional model points where there 
are large errors. Points with the smallest coefficients 
are eliminated from the set. A new model is computed and 
the process is repeated. This procedure can be repeated 
until the desired accuracy is reached or until additional 
iterations fail to produce increased accuracy. This method 
is considerably better than selecting regular spaced points 
and densifying them for more accuracy or selecting the maxi- 
mum number of points at one time by examining the bottom 
features. The iterative method of model point selection is 
adaptable to automation since there is relatively little 
subjectivity involved. 
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Adjacent survey models or partial survey models can 
be joined to produce a continuous and imambiguous repre- 
sentation of the bottom in the junction area. The best 
method is to use a v^eighted sum of the model values in an 
overlapping area. The Hermite pol 5 niomial function should 
be used to produce the weights. If the two models contain 
common points in the overlapping area, a good junction can 
be made without computing a v/eighted sum. This method 
could produce a slight discontinuity along the boundary 
but it was not apparent in the test runs for this project. 
The multiquadric technique with iterative selection 
of model points and Hermite pol 3 momial jtuictions produced 
good models of foxir data sets. Approximately 20 % of the 
survey data points were reqtiired for a model of A^lke Bay, 
Alaska, v/here there is tremendous bottom irregularity. 

Only 10% of the data points were required for a model of 
the Gulf Coast v/here the bottom shov/s little variation. 
Approximately 15% and 17% were required for the Horro Bay 
and Monterey Bay models where there is more irregularity 
and moderately sloping bottoms. 
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